in 

a^ 

Cellular Automaton for Realistic Modelling of 

r^ Landslides 

(N 

! E. Segre 

^ ' Istituto di Cosmo-Geofisica del C.N.R., Torino, Italy 

(N ' 

O ' C. Deangeli 

^ . Politecnico, Dipartiinento di Georisorse e Territorio, Torino, Italy 

O 

a 

I ■ Camera-ready Copy for 

r^ , Nonlinear Processes in Geophysics 

O '. Manuscript-No. 940030, version 2 

> 

k> , Offset requests to: 

H I Enrico Segre 

' 10/1994 - 3/1995: RIMS, Kyoto University, 606 Kyoto, Japan; 

4/1995 - 3/1996: Dept. of Civil Engineering, Technion, Haifa, Israel 
e-mail: segre@to.infn.it, segreOkurims.kyoto-u.ac.jp 



Journal: Nonlinear Processes in Geophysics 
MS No.: 940030, version 2 
First author: Segre 



Abstract. A numerical model is developed for the sim- 
ulation of debris flow in landslides over a complex three 
dimensional topography. The model is based on a lat- 
tice, in which debris can be transferred among nearest 
neighbors according to established empirical relation- 
ships for granular flows. The model is then validated by 
comparing a simulation with reported field data. Our 
model is in fact a realistic elaboration of simpler "sand- 
pile automata" , which have in recent years been studied 
as supposedly paradigmatic of "self-organized critical- 
ity" . Statistics and scaling properties of the simulation 
are examined, and show that the model has an intermit- 
tent behavior. 



1 Background 

Landslides are among the processes responsible for the 
evolution of the shape of the earths surface (Scheideg- 
ger, 1991). They constitute a major natural hazard, and 
therefore their prediction and modelling are very impor- 
tant in the environmental sciences. 

The approach generally undertaken in soil engineering 
for the study of landslides concentrates on the phase of 
the slope failure, i.e., the analysis of the conditions for 
the breakdown of an elastic soil (see, e.g., Hutchinson, 
1988). This approach, despite being fundamental for the 
processes referred in geotechnics as falls, topples, slides 
or spreads, which actually initiate a landslide, does not 
apply after the complete liquefaction of the soil mass, 
and does not describe the subsequent debris flow. 

The establishment of constitutive equations for debris 
flow is, however, not trivial. Debris flow is in general 
described as a heterogeneous suspension of granular soil 
particles in water, in proportions and with grain size 
distributions that may largely vary. The overall motion 
is clearly too complex to be described by tracking each 
particle, and therefore the mixture is viewed as a con- 
tinuum. The granular nature of the debris gives rise to 



a non-Newtonian fluid-like behavior. Since the seminal 
work of Bagnold (1954), a huge number of rheological 
models have been proposed in the literature for differ- 
ent possible regimes of these flows (see e.g. Chen, 1987; 
Sauret, 1987; Takahashi, 1991). Both the nomenclature 
and the models vary largely from author to author, and 
are often differentiated according to the initiating causes 
or the kind of soil involved in the slide (Pierson and 
Costa, 1987). 

Since landslides have catastrophic effects, and cannot 
often be studied effectively under controlled conditions, 
one is led toward their computer simulation. To the 
best of our knowledge, computational debris flow models 
previously given in the literature have been quite sparse 
(Sassa, 1988; VuUiet and Hutter, 1988; Takahashi, 1991), 
and were based on finite elements implementations of 
the assumed fluid equations. 

In contrast, in this paper we build a model which in- 
cludes a stopping and re-initiation mechanism that can 
give rise to criticality Our aim, apart from that of to 
simulate effectively landslides, is to reconnect with re- 
ported results for the so-called cellular automata. 

Cellular automata are computational lattice models in 
which each reticular site can be in one of a limited num- 
ber of states, whose evolution is solely determined by 
local rules. In recent years these models have been em- 
ployed in many different physical computations. Some 
of them, later called lattice gases, have been proved 
to conform in the thermodynamic limit to important 
macroscopic differential equations (Frisch et al., 1986; 
Chen et al., 1992, 1993). In this respect cellular au- 
tomata can be translated into algorithms which are sig- 
nificantly more efficient than corresponding finite dif- 
ferences integration schemes, and are straightforwardly 
applied to complex boundary geometries. 

Moreover, a particular cellular automaton, namely 
the "sandpile model" of Bak et al. (1987), has been in- 
troduced in connection with avalanche dynamics. That 
automaton has been used as a toy model to explain 
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the occurrence of power-law statistics in natural pro- 
cesses, via a supposedly universal feature which has been 
named self-organized criticality. Since the original pa- 
per, there has been much discussion in the literature 
about the real significance of that model, and about 
its applicability to real granular flows. Various elabo- 
rations of the original automaton of Bak et al. (1987) 
have been reported since then. Some of these elabora- 
tions aim to achieve more realistic effects, and gener- 
ate different statistical and scaling properties. Recent 
contributions to the subject include for instance those 
by Dhar (1990); Prado and Olami (1992); Ding et al. 
(1992); Christensen et al. (1992); Krug (1993); Theiler 
(1993); Rosendahl, et al. (1993); Carlson, et al. (1993); 
Morales-Gamboa et al. (1993); Chhabra et al. (1993); 
Maslov and Olami (1993); Frette (1993). Experiments 
with real sand grains, as well with other granular mate- 
rials have been carried out for example by Nagel (1992); 
Evesque et al. (1993); Nori and Pla (1993), and com- 
pared with the expectations from theoretical sandpile 
automata. A comprehensive review of the experimental 
and computational investigation on granular -and specif- 
ically sand- flows can be found in Metha and Barker 
(1994). A descriptive contribution by Noever (1993), 
which considers statistics of landslides in a mountain 
region, is also to be mentioned in this context. 

The aim of this paper is to propose an effective cel- 
lular model for debris flow, combining the schematicity 
of sandpile models and more appropriate geotechnical 
empirical laws, and to look at its statistical behavior. 
To our knowledge, a similar approach has not yet been 
undertaken for debris flow in any other place than the 
note of Digregorio et al. (1994), which appeared after 
the first version of this paper was submitted. However, 
a much simpler cellular automaton is considered there, 
and the attention is mainly posed on its parallelization. 

The paper is organized as follows: the rules of our au- 
tomaton are derived from established debris flow laws 
in Sect. 2; the scheme is numerically tested in Sect. 3; 
a real event is simulated in order to calibrate the model 
in Sect. 4, and finally a multifractal analysis of the re- 
sulting time series is reported in Sect. 5. 



2 Computational model 

2.1 Lattice setup 

We are primarily interested in describing the fiow of a 
heterogeneous and scarcely coherent debris, such as that 
of a stony discharge, which takes part over a rigid and 
non erodible substratum. We furthermore assume fully 
mixed debris {mature debris flow in the nomenclature 
of Takahashi, 1991), i.e., a homogeneously saturated de- 
bris layer. The appearance of surface fiowing water, the 
shear stability of different strata, etc., will not be upon 
question here. Our model docs not account for varia- 



Table 1. Data values interpolated for the constitutive relations 



(silt) 


S2 
(gravel) 


S3 
(boulders) 


ps 

g/cm'^ 


degrees 


100% 


0% 


0% 


1.8 


12° 


0% 


100% 


0% 


2.0 


30° 


0% 


0% 


100% 


2.2 


35° 



tions of the debris properties in the vertical direction z, 
i.e., we adopt a vertically averaged description. This al- 
lows us to schematize to some extent the real dynamics. 
The landslide site is discretized in elementary cells of 
flnitc size. In each one the state of the system is specified 
by the values of some representative quantities. These 
include: the height of the impermeable rigid substratum 
r; the nonnegative amounts of water w, of gas (or void 
fraction) g and of granular solids s in the cell. In order to 
account for debris properties that depend on the gran- 
ulometry, the solid content is classified by size. Three 
size classes have been presently employed, representing 
silt-clay, sand-gravel and boulders. All the contents are 
given as partial heights (volumes/base area of the cell), 
so that the top height in the cell i is given by 



z{i) = r{i) + w{i) + g{i) + ^ Sk{i) , 



(1) 



fc=i 



and the volume of its content is V = S[z{i)—r{i)\, where 
S is the base area of the cell. 

The density of the solid part ps and the "friction an- 
gle" (f) depend on the composition of the debris. In first 
approximation they are calculated by linear interpola- 
tion of the data points given in Table 1, as 



Ps{i) 



T,l=iPkSk{i) 
Efe=iSfe(«) 



Efc=iSfc(«) 



(2) 



i.e., the two quantities are expressed as weighted means 
of the ones of the individual solid components. 

More sophisticated empirical laws, detailing different 
soil classes, could be employed in a future refinement; 
the form of ps used here allows to impose the conser- 
vation of mass simply by conserving the volume of the 
transported material, and fits satisfactorily to reported 
soil data. 

The embedding fluid is assumed to be clear water, 
with density pf — Ig/cm^; the weight of gas is ne- 
glected. The solid concentration C and the density p 
of the mixture in cell i arc thus 



Ciz) 



J2l=iSk{i) 



g{i)+w(i) + J2k=iSk{i) 



p{2)^Ciz)ps{i) + {l-C{t))pf. 



(3) 



(4) 



Each cell is considered connected with a number b of 
its nearest neighbors. At each discrete time step some 
material can be transferred to those neighbors, accord- 
ing to rules which take into account the local slope in the 
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Fig. 1. Diflfercnt lattice geometries based on a cartesian mesh: 
a) square, b) hexagonal, c) octagonal, d) triangular, e) 16-stellar. 
Note that in case b) the white points constitute a shifted inde- 
pendent lattice. 



neighbors directions. The whole lattice is thus viewed 
as a network of elementary slopelets, each one with its 
characteristic slope due to the differential topography 
and debris accumulation. 

Different lattice geometries, corresponding to different 
choices of meshes and sets of b first neighbors, can be de- 
vised. The most obvious arrangements are the Cartesian 
square lattice {b = 4) and the hexagonal lattice (6 = 6). 
Other staggered stencils are easily implemented (Fig. 1). 
A geometrical realization of the top surface may have to 
be dropped, because elementary slopelets cross them- 
selves as the stencil is shifted. Such arrangements may 
be, however, effective in smoothing potential discretiza- 
tion instabilities (see below in Sect. 3); the set of local 
slopes is then interpreted as a pure property assigned to 
the lattice. In such cases the fact that cell j is a neigh- 
bor of cell i may not even imply the converse; consider 
for instance what is obtained by shifting the pattern of 
Fig. Id along the lattice. 

We have to resort to a two-neighbors partition rule for 
the transfers, in order to achieve lattice isotropy (Fig. 
2). By this we mean that the flow along an arbitrarily 



Fig. 2. Schematic of the debris transport process (depicted for 
square cells, but easily transposed to other geometries). 



inclined plane must not depend on the orientation of 
the lattice. To this extent a local slope angle 9{i,v) is 
defined in the cell i for each pair of its adjacent neighbors 
{jv,jv+i} {v = 1, . . . , 6 modulo b). The steepest descent 
vector p(i,v) on the plane passing for the three center 
points is 



Pihv) = {jy,Az^{i)) A (j;+i,A2;„+i(i)), 



(5) 



where fy denotes the two-dimensional lattice vector ori- 
ented from node i to jy and Az„(i) — z{jy) — z{i). The 
slope angle is then 



(i, v) = tan" 



-1 \Pxy{hv)\ 
-p{i,v) -z 



(6) 



where zis the unit versor in the z direction, and Pxyi^i, v) 
the projection of p(z, v) onto the xy plane. The outflow 
rate q{i, v) from the cell i toward the sector {jv, jv+i} is 
assumed to be parallel to p{i, v) and function of 9{i, v). 
For this purpose Pxy{i, v) is decomposed as 



Pxy{i,v) ^pi{i,v)ey + p2{i,v)ey+i , 



Pi{i,v) 
P2{i,v) 



1 - {cy ■ Sy+iy 



\p{i,v)\ 
\p(i,v)\ 



1 



{Cy ■ Cy + l)^ 



(7) 



(8) 



where ey = fy/\fy\ is the w-th lattice versor. 

We impose volume and, consequently, mass conserva- 
tion separately for each solid class and for water. The 
gas balance will be taken into account differently, since 
free exchange can take part with the atmosphere. A 
more elaborate model could eventually consider either 
fragmentation processes, in which some amount of the 
larger granular classes may break and add up to the 
smaller ones. Alternatively, there may be more overlaid 
strata, each one with its composition, with possible mu- 
tual exchanges. In any event, we obtain a model which 
has, in the sandpile terminology, local and limited rules. 

Energy and momentum conservation are not enforced 
here. This is consistent with modelling a process which 
is dissipative in character, both because of the increas- 
ing dispersion and liquefaction and because of viscosity. 
The assumption that the debris previously in motion 
can suddenly be stopped within the space of a single 
cell, depending only on instantaneous local conditions 
(as modelled below), is equivalent to the assumption 
that kinetic energy is readily dissipated, gravity being 
the main energy source. 

Initial conditions for the model are imposed specify- 
ing the site topography and the debris distribution at 
the initial time. Since our primary aim is to model de- 
bris discharges, which possess typical time scales rang- 
ing from seconds to minutes, we will consider the total 
water amount as flxed. Even if the primary cause which 
triggers landslides is usually an extreme rainfall, the wa- 
ter content of the soil varies on a longer time scale, and 
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need not to be externally varied during the simulation. 
For the same reason, we think that permeation through 
the soil has no time to occur in this regime of flow, and 
we do not consider the use of conductivity laws, such as 
Darcy's law. Net in- or out-fluxes in particular source 
or pit cells can, however, readily be included into the 
model, and could be required as boundary conditions. 
This would be the case, for example, if the computa- 
tional domain had to be truncated up or downstream, 
or if permeation through the substratum intervened sig- 
nificantly into the water balance. 

Boundary conditions are easily implemented on any 
set of cells: to realize an open boundary it is just suf- 
ficient to force the content of these cells to be always 
null, while to achieve closed boundaries one simply sets 
to zero any resulting out-flux. Neighbors relative posi- 
tions could be calculated modulo the lattice dimensions 
to achieve periodic boundary conditions, if ever neces- 
sary. 

2.2 Deriving automaton rules from empirical flow laws 

Evolution rules for the automaton are imposed consis- 
tently with established empirical laws for debris flow. 
These rules must account for the initiation ("toppling" 
in the sandpile literature), motion and stopping of the 
debris mixture. As a general issue (see Chen, 1987), 
we note that often the rheological models underlying 
the empirical laws may not be completely correct. The 
model equations may not be able to account for the ob- 
served values when microscopic quantities (such as the 
actual grain sizes) are substituted. The proposed for- 
mulas can nevertheless sometimes fit real data, if one 
relaxes the assumptions that were made to derive them 
or introduces ad hoc empirical factors. This fact encour- 
ages us to make use of rather empirical rules, regarding 
the approximate agreement with real data as the main 
validation. 

Under the conditions exposed at the beginning of the 
previous paragraph, the debris flow is in the so called 
granulo-inertial range (Seminara and Tubino, 1990; Ta- 
kahashi, 1991). In this dynamical range the forces due 
to the turbulent motion of the embedding fluid are neg- 
ligible with respect to the coUisional stresses of the solid 
fraction. The macroscopic behavior of the granular mix- 
ture is assimilable to that of a dilatant viscous fluid, i.e., 
a non-Newtonian fluid with a quadratic dependence of 
the viscous forces upon the shear stresses. 

We model the elementary slip processes solely with a 
state of rest and a state of flow. A rationale for doing 
this is that the distance required for complete liquefac- 
tion of a debris mass is short (see Takahashi, 1991, Fig- 
ure 3.6). Therefore, in a simplified model, we neglect the 
details of the (slope failure) initiation phase of the move- 
ment. We will also make use of the assumption, perhaps 
strong, but adequate for our purposes, that the elemen- 
tary fiows along the lattice take the form of granulo- 



inertial stationary debris flow in uniform rigid bottom 
channels. 

2.2.1 Initiation 

A condition for the occurrence of debris flow on a slope 
plain is given by Chen (1987), Eq. (42): steady equilib- 
rium flow takes part if the slope angle 9 satisfies 



(9) 



^^^sin0<tane<-^^^^. 

In Eq. (9), (p is called the "static friction angle" and 
fii and fj,2 represent respectively a "consistency" and 
a "cross consistency" coefficient of the squared shears 
in the expression of the dilatant viscous stresses. Vari- 
ous forms for the dependence of these coefficients on the 
concentration and on the granulometry are proposed in 
the literature, leading to different results and leaving 
much space to ad hoc fits. For example these coeffi- 
cients are estimated by Chen (1987) to depend on the 
concentration as ^(1 — C/C*)"^-^*^', as is suggested by 
Einstein's equations for granular viscosity, with C* equal 
to the maximum admissible concentration and A a con- 
stant. By referring to the original Bagnold's work (Bag- 
nold, 1954), Takahashi (1991) instead expresses /ii = 
aps(\dY sin^d, with a an empirical constant, d the solid 
particle diameter, A the linear concentration and (pd the 
"Bagnold quasistatic stress angle" . 

Some disagreement also exists in the literature about 
the use in Eq. (9) of <j) rather than c^^, which are found 
experimentally to differ by a few degrees. The lower 
bound on 9 coincides with the Bagnold estimate if it is 
assumed that tant/)^ = sin0 in the quasistatic limit. If, 
furthermore, ^2 — —cipsi^d)^ cos(j)d is assumed, Eq. (9) 
would tell that equilibrium steady uniform flow occurs 
only for 9 — (pd] that would however lead to other in- 
consistencies, as pointed out both by Chen (1987) and 
Takahashi (1991). Since in our study the angle 4> is em- 
ployed rather heuristically, we will not need to be more 
precise. More care would be needed if the appropriate 
value was to be compared with specific in situ measure- 
ments. 

In our lattice model we neglect the upper bound and 
simply translate Eq. (9) into the rule: the flow initiates 
on every lattice sector (i,u) for which 



iSi\i9(i,v) > ^-^^--^tan( 
PW 



(10) 



The intervening p{i), 9{i,v) and (j){i) are calculated in 
each cell as explained in Sect. 2.1. For every cell of the 
lattice, all the branches to neighbors j„ such that Eq. 
(10) is satisfied will be called critical directions, and flow 
is assumed to move material to those neighbors only. 

2.2.2 Flow 

In the stationary flow of a granulo-inertial debris layer 
over a rigid bottom, the flow velocity u has a vertical 
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dependence on z like (Chen, 1987, Eqs. 46-49): 



2 p^gsinO 
uiz) — —\ 



3/2 



(^0 - z) 



3/2 



where g is the gravity acceleration, p* = 7p is a cor- 
rected density, with < 7 < 1. In Eq. (11), zq is the 
height from the substratum at which the flow velocity 
becomes independent from z, e.g., where the shear stress 
equals the vertical component of the weight of the up- 
per layer. For flows of incohesive granules, zq = h (the 
whole layer of debris moves according to the law) and 
7^1. Vertical averaging of Eq. (11) gives then the 
cross sectional velocity U of the entire layer (Takahashi, 
1991, Eq. 2.3.4): 



2 / gph'^ sin ( 



(12) 



This formula is consistent with dimensional estimates 
based on the classical Manning or Chezy forms of the 
stresses (Trowbridge, 1987; Scheidegger, 1991), which 
account for the presence of the square root of psin^. 

If the concentration of the solids or the granulomet- 
ric distribution are allowed to vary within the debris 
depth, as it is frequently observed, then the descrip- 
tion becomes more complicate. A number of peculiar 
phenomena are observed in non-homogeneous granular 
flows; for instance the larger boulders tend to float over 
the bulk of the debris, can be transported farther apart 
and accumulate along the fronts. Qualitative explana- 
tions for the observed levitation of the larger grains are 
found in Savage and Lun (1988) and Takahashi (1991), 
§4.4. 

As for the concentration of the flowing debris, a hint 
could be taken from the equilibrium concentration for 
the steady uniform inertial debris flow, which is 



C = 



Pf 



tan 6' 



Ps — Pf tan (j) — tan 9 



(13) 



(Chen, 1987, Eq. 62; Takahashi, 1991, Eq. 2.3.9, though 
inconsistent for very steep channels). The use of such an 
equation would imply to assume that exactly as much 
water takes part in the motion as is required to maintain 
the steady uniform debris flow. Since it is not assured 
that the required amount of water is actually available 
in the cell, we prefer to treat the different components 
as if they were moving independently. At the same time 
we can mimic the faster velocity of the larger grains. 

We therefore implement a variant of Eq. (12) on the 
lattice. We let quantities q of material flow out of the cell 
i toward its critical neighbors and evaluate q by putting 
z{i) — r{i) as h and 0{i,v) as 6 in Eq. (12). We assume 
that the concentration and the density of the flow are 
those of the source cell i. The slope 0{i,v) would also 
decrease during the elementary motion. We would have 



to account for a temporal dependence of h(t), obtained 
by solving the differential equation 



(11) dh{t) _q{C,t) 



lgp{t)h{t)'^sme{i,v,t) 



dt 



S 



MiW 



(14) 



in which P is an appropriate effective width (hydraulic 
radius), and also the concentration C(t), the density 
p{t) and the slope inclination depend on time. In order 
to keep the model as simple as possible, we assume that 
the q{C, Az) are constant during the time step, and rely 
on the first order approximation for small values of the 
time step Ai. 

Moreover, we compute independent q^ for each gran- 
ular fraction (fc = 1,2,3) and for water (fe = 0), and 
assume complete mixing before and after the transfer. 
The simplified evaluation is thus 



qk 



{i,v) = Q k At y^ ph^ sin 9 {i 



(15) 



with Qk constants, conglobating and approximating all 
the neglected factors of Eq. (14). The ratios of the Qk 
entirely account for the faster propagation of the larger 
grains. Water is treated as a granular class, since we 
put w{i) ~ so(J) and use Eq. (15) with a rate constant 
Qo is of the same order of the Qk- However, water is 
differentiated from the solids since it plays a different 
role in the evaluation of p, of and in Eq. (10). 

As mentioned before, the outflow from each sector 
{i,v) of the cell reparts among two adjacent neighbors 
proportionally to pi andp2 computed from Eq. (8). The 
time advancement rule will be 



Sfc(i) -^ Sk{i) -[pi{i,v) +p2{i,v)]qk{i,v) 
Skijv) -^ Skijv) +Pi{i,v)qk(i,v) 
Sk{jv+i) -> Sk{jv+i)+P2{i,v)qk{i,v), 



(16) 



where it is intended that qk{i,v) = if the w-th sector 
is not critical. Upper limits on qk{i, v) prevent the cells 
to be more than depleted within the time step: after all 
qk{i,v) have been calculated for the cell «, 



qk{i,v) -^ qk{i, v)nun 1 



Sk{i,v) 



(17) 



and only then the Sk{i) are updated. The rule Eq. (16) 
with the correction of Eq. (17) is iterated over all i and 
all V. All the computed qk{i,v) are stored before the 
simultaneous updating of the lattice. In this way the 
computation may also be completely parallelized. 

We treat the void fraction in a similar way. Void, 
which has to be taken into account for the total volume, 
will in our simplification be class —1, moving around 
the lattice according to Eq. (15) with its own rate fac- 
tor Q-i. However, void is not conserved: experience 
shows that shaken granular material tends to unpack 
and hence to augment its volume (see e.g. Metha and 
Barker, 1994). To keep the model schematic, we let the 
running material increase proportionally its void ratio. 
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Table 2. Maximum allowed void ratio for each class 





Si 


S2 


S3 


1 f^max 


1.5 


0.835 


0.3 



till a maximum value Cmax- As a consequence, the void 
is propagated to the neighbors with the rule 



S-l{jv) 
S-lijv+l) 



s-i{i) - [Pi +V2]q-i{i,v) 

3 

q-i +a'^qk 



s-iijv) +Pi 

S-l{jv+l) +P2 



k=0 



5-1 + '^^qk 

k=0 . 



(18) 
ii,v) 

(*,f), 



where a is a parameter to be calibrated, and the qk{h v) 
have been calculated previously with Eqs. (15)-(17). 
Subsequently the upper limitation is enforced: 



iW 



s-i(i) 



k=l / 



(19) 



The limitation involves the maximum void ratio emax^S) 
admitted in the destination cell after the transfer. This 
ratio is calculated as that of a composite soil of corre- 
sponding composition, that is, by weighted average of 
the values given in Table 2 (Lancellotta, 1987, p. 7). 
The maximum void ratio is computed according to the 
solid fraction only. 

The time advancement proceeds by discrete steps. 
The physical value of Ai will be a consequence of the 
dimensional change that is allowed to take part in the 
system during a single step. This is in turn determined 
by the Qk^ which should be small enough to ensure the 
near constancy of h(i) within a time step. An estimate 
of the temporal evolution of the outflow from a single 
cell is obtained by numerical integration of Eq. (14) and 
(15) as At -^ dt. In such case the numerical integration 
shows an asymptotic relaxation of /i(i) oc t~^/^. 

Instantaneous flow rates q(i) are evaluated in each cell 
by vector summing all the incoming flows, that is. 



'TIO = E 



j(, 



'z-i(.?tO + (i + a)I]*0'«' 



fe=0 



where j^ labels all the cells for which i is a first neigh- 
bor. Local instantaneous velocities are derived from 
these flow rates by dividing for At and appropriate ge- 
ometrical width factors. 

No further assumptions are made about the dynam- 
ics of the elementary movement. The cell forgets the 
details, and smoothes the dishomogeneities after the ad- 
vancement step has been made. The elementary transfer 
is not resolved below the Ai scale. 

We note that despite our scheme does not directly 
account for forces, it is not completely non dynamical, 
since it employs flow laws that have been derived with an 



underlying rheological model, and it still retains features 
of the viscous flow. 

2.2.3 Stoppage or continuation 

Very few quantitative arguments are given in the litera- 
ture for the exact calculation of the stoppage conditions 
of a debris flow which runs on an undercritical slope. 
At least to our knowledge, the exact determination of 
the stoppage distance of a front is still an open problem 
in the general case (Takahashi, 1991, § 5.1). What it 
is at least established (Jaeger and Nagel, 1992; Metha 
and Barker, 1994) is that the repose angle at which col- 
lapsing dry granular material settles is by some degrees 
lower than the critical initiation angle. 

We included in our model a simple local rule that 
accounts for the dynamical lowering of that angle. For 
this, a one-step memory is added to the model. If the 
cell received material at the previous time step, then still 
enough kinetic energy may be available to sustain the 
flow in a non critical direction. In this way a fast enough 
flow can even climb counter-slopes, as it is sometimes 
observed. The slope direction vector p{i,v) appearing 
in Eq. (6) is dynamically corrected in 



p{i, v) -^ p{i, v) + l3q{i, t - At) , 



(21) 



where q{i,t — At) is the volume inflow in cell i at the 
preceding time step, and /3 is another parameter of the 
model. The flow then continues if the corrected 9{i,v) 
satisfies Eq. (10). In other words, it is assumed that 
an "effective slope" is seen by the running debris. The 
argument could be related to energy line considerations 
(see e.g. Sassa, 1988). A more elaborate model could 
take a different threshold into account, or allow for a 
partial deposition of material. 

2.3 Summary of the computational model 

For the algorithmic implementation of Eqs. (16), (17), 
(18) and (20), it results convenient to organize the vari- 
able information into two arrays. The first one is Sk{i) 
(—1 < fc < 3), and the second is the transfer array 



(20) Gk{i)^sk{i,t + At)-Sk{i,t). 



(22) 



in which the differences result from the above mentioned 
equations. In other words, it is more convenient to track 
the in-fluxes coming from all j(, toward the cell i, than 
to track the out-fluxes going to each j„. 

The scheme requires, at each time step, the following 
operations to be performed: 

— the state of each cell i is defined by r{i), Sk{i), and 
by q(i) at the preceding step; 

— z{i), p{i) and 0(i) are computed by Eqs. (1), (2-4); 

— 0{i,v) is computed by Eqs. (6) and (21) for all 
sectors of i; 
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— the critical directions are found by Eq. (10); 

— the qk{i,v) are calculated by to Eq. (15-16); 

— thresholds are enforced according to Eq. (17); 

— the void increase is computed by Eq. (18); 

— all gains and losses are included in Gk{i)', 

— the boundary cells conditions are enforced on Gk{i)', 

— the influx rates q(i) are computed by Eq. (20); 

— the whole lattice is updated adding Gk{i) to Sk{i)', 

— the void limitation, Eq. (19), is enforced. 

3 Testing the numerical model 

3.1 Numerical stability 

The implementation of the automaton with continuous 
state variables renders the explicit time advancement 
scheme only conditionally stable. We have observed, 
for example, the appearance of an oscillatory instability 
pattern. This can occur whenever, at a given time step, 
an exceeding quantity of material is transferred from a 
cell to its neighbor(s). If the material is then not drained 
away rapidly enough, then at the following time step the 
slope between these cells may result reversed, so that 
the the first outflow just returns into the original cell. If 
the mechanism is allowed to repeat, the back and forth 
transfer of matter gives rise to persistent oscillations, in 
which odd and even cells are alternatively depleted and 
flUed. Averages either over a few time-steps or over a 
few cells remain well behaved in presence of this instabil- 
ity; but the unphysical oscillatory background velocities 
may amplify if [3 is large enough. With our set of rules, 
enforcing ulterior thresholds to prevent this instability 
would be cumbersome. The best cure to this problem 
is a refinement of the time step (i.e., a rescaling of all 
Qk), which scales down the elementary outflows. 

The actual criterion for stability can be outlined in 
the following way: consider cell 1 and 2 upon a flat 
bottom, spatially separated by I, which contain at time 
t respectively hi and /12 heights of equal debris. Let the 
resulting slope angle be such that 



It must result tan^j 



> — tan (j) in order to prevent a 
return of debris from 2 to 1 , and therefore 



hi — h2 
tant'i = > tan( 

A partial height 



dAt - 0(1 + a) At J phi sin 01 



(23) 



(24) 



is then transferred from 1 to 2, where Q is an average of 
the Qk- Thus at time t+At the effective slope, including 
the advection contribution, becomes 



At < 



/(tan( 



tan( 



2d 



21 



(26) 



tan 62 = tan 61 



2dAt (3d 



(25) 



An upper bound on At in terms of typical ft,, Z, Q, p can 
thus be derived. Since the outflow from one cell may 
be shared among more neighbors, this bound may be 
lowered by a factor up to h. In this sense lattices with 
higher coordination number may be numerically more 
stable. 

3.2 Isotropy 

In order to claim that our model faithfully reproduces 
the physical process, we must ensure that it has an 
isotropic behavior, at least in an averaged sense. This 
implies that the flow computation must result as inde- 
pendent as possible from the choice of the lattice geom- 
etry and from the orientation of the lattice with respect 
to the landscape features. 

The two-neighbors transfer rule introduced in Sect. 
2.1 guarantees a local spatial isotropy. We note that 
it is the use of flow laws such as those described in 
the previous sections, that renders necessary the two- 
branches mechanism. It can be easily checked that if a 
simpler single branch rule was employed, then the model 
would have been isotropic for any lattice geometry only 
if g oc Az, and no threshold on the minimal slope angle 
was enforced. 

A theoretical verification of the isotropy requirement 
is not easy, because the flow evolution is determined 
by the continuous buildup of the landscape. It is not 
easy, for instance, to verify whether the flow reaches at 
the same time equally distant cells, regardless of their 
positions on the lattice. A detailed calculation of the 
transport between two sites would require to consider 
the contributions through every possible path connect- 
ing them. None of these paths is independent, because 
the material that accumulates along the paths alters 
the global landscape. The evaluation of the expecta- 
tion or of the correlation values of the flux in different 
points or directions would be similarly affected. The ar- 
guments normally quoted for the isotropic behavior of 
lattice gases onto the hexagonal lattice (referring, e.g., 
to Frisch et al., 1986) do not apply to the present case, 
for the much too different setting of the problem. 

Some isotropy in a time averaged sense is however 
numerically seen to be recovered. We choose to verify 
the isotropy of the model by considering the transport of 
matter in arbitrary directions. A number of numerical 
experiments can be carried out to this extent. We report 
two of them that convince us of the isotropic behavior. 

To begin, it is verified that the flow rate is indepen- 
dent from the slope line orientation. Figure 3a shows 
the flow vectors obtained in a test simulation of the dis- 
charge of an uniform debris layer over a conical pedestal. 
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Fig. 5. Roll wave pattern in a numerical simulation of a channel 
debris flow with constant upstream alimentation. Grid resolution 
is 61 X 61. 



Fig. 3. flow vectors of an initially uniform distribution over a 
conical pedestal, after the first time step; b) scatter-plot of the 
velocity vectors of case (a) (the few off-circumference points rep- 
resent either boundary cell rates, or are discretization artifacts) 



The resulting flow is expectedly radially symmetric, and 
the flow vectors show, apart from effects due to the dis- 
cretization, a radial distribution, as can be seen from 
the scatter-plot in Fig. 3b. 

Another test, inspired by Prigozhin (1993), is to check 
that material emitted by a single source cell over a flat 
bottom accumulates in a circular based conoid. This 
is also achieved, and shown in Fig. 4. The parameters 
used in these two tests are the same used in the case 
simulation described below in Sect. 4. 

3.3 Roll waves 



Channelized steady debris flow can develop an undula- 
tory free surface instability, namely the appearance of 
roll waves, which are often observed in real events. Con- 
ditions for the development of such instability, which im- 
ply an evaluation of the shear and of the bottom stresses 
of the debris, have been examined by Trowbridge (1987) 
for the Newton and for the Bingham fluid model. We 
were able to achieve a similar-looking pattern in the 
numerical simulation of a constant slope channel, with 
upstream constant alimentation and downstream open 
boundary conditions. The debris layer thickness after a 
long simulation time is presented in Fig. 5; here again 
the same parameters of the case studied in Sect. 4 have 
been used. 



Fig. 4. Accumulation conoid produced by a point source on a 
flat topography, after 40000 time steps of 0.02 sec with /3 = 0. 



4 Case study of a real landslide: validation and 
back analysis of the model 

The agreement of the computational result with a real, 
full size landslide validates our model. Since during the 
time of this work we had no direct access to a fully mon- 
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itored event, we had to resort to a case history found 
in the literature (Jiang Yun et al., 1991). This event 
was chosen among others because of the extended in- 
terested region and the relative variety of phenomena, 
which included also counter-slope flows. We were able 
to reproduce the major features of this event from the 
data presented in that report, even conjecturing some 
of the missing information. 

The landslide detached from the Xi Kou Mountain in 
the Sichuan province of China, on July 10, 1989. The 
slope where the landslide took place was composed by 
very weathered and fractured strata of limestone, mar- 
lite, dolomite and shale. The bedding had a dip opposite 
to that of the slope. The ridge of the slope was 1012?ti 
high, while the Xi Kou town is located on a depression 
at a height of 298?7i. Between the heights of SOOtti and 
600™ the slope was initially covered by a debris blanket 
of triangular shape. This debris blanket had a thickness 
varying from a few centimeters at the top of the triangle 
to ten meters at its basis. On its lower side there were 
deposits constituted by finer materials, which probably 
came from previous landslides. The blanket lay on a 
natural platform, the Maan Pin platform, constituted 
by shale and dolomite strata. That platform was about 
200m long, and it was covered by a thick 10?n clay layer. 
The platform was weakly dipped downwards; at its end 
began a steeper slope, with a difference of level of about 
40m. After this jump there was a long channel inclined 
on average of 19°. On the right side of the platform, a 
hill of finer material was also interested by the develop- 
ing flow. The reported landslide was probably triggered 
by the heavy rainstorm occurred on the preceding days. 
The water, flowing into the rock fractures, saturated the 
rock mass and therefore caused an overpressure under 
the clay layer. The overpressure decreased the effective 
stresses and led to the slide movement. 

It is possible to divide the landslide in two phases. 
The flrst one is the slide of the clay layer and of the 
debris blanket which lay on the platform. After running 
out of the Maan Pin, this material fell into the 40m deep 
cross channel and decomposed into debris flow. The sec- 
ond phase of the landslide is a proper debris flow, which 
split in two branches. The main one ran downward along 
the gully and accumulated in a fan at the slope toe. The 
mass of the other branch eroded the hill, ran along a 
right gully and accumulated at the gully curve. 

The total volume of the material involved in the land- 
slide was estimated to be 1.000.000 m^, about ten times 
larger than that of the initially detaching mass. This 
fact indicates that, during its flow, the debris both in- 
corporated other material along the path and increased 
its voids ratio, thus increasing its total volume. 

Our numerical simulation treated the entire process as 
a debris flow. This has been necessary in order to over- 
come the difficulty of establishing accurately the shape, 
the quantity and the velocity of the material at the pass- 
ing from slide to debris flow. 



Fig. 6. Topography of the Xi Kou landslide, as obtained from 
Jiang Yun et al. (1991). Bold lines represent the isolines acquired 
from the original map; the reconstructed isolines and the contours 
of the reported event zones are shown in the upper projection. 



The topography of the site was reconstructed by ac- 
quiring the planimetry data reported by Jiang Yun et 
al. (1991). A total of 546 coordinate points on the map 
was tracked; an interpolating surface (Fig. 6) was then 
computed in order to obtain the topographical surface 
heights at regularly spaced lattice positions. An inter- 
polation library routine based on a radial functions fit 
(PV-WAVE, 1993, routine RADBF) was employed. The 
resulting lattice was a 61 x 151 rectangular grid of 10m 
spaced points. The thickness and the composition of 
the debris layer were defined in every lattice cell accord- 
ing to the description of the site given by Jiang Yun 
et al. (1991). Three zones of the site were filled with 
mobilizable material, with the different compositions re- 
ported in Table 3. The water content was determined 
by assuming that, initially, the material was completely 
saturated: its entire void volume, which is found multi- 
plying the solid volume by cq, was assumed to be filled 
by water. The "in situ void ratio" cq < Cmax was deter- 
mined from empiric tables (Lancellotta, 1987). 

The results of the numerical simulation are depicted 



Table 3. Initial distribution of mobilizable material for the Xi 
Kou slide simulation 





Triangle 


Platform 


Hill 


Shape 


pyramidal 


10m layer 


14m conoid 


Sl % 


33.1 


35.7 


24.3 


S2 % 


49.7 


23.8 


36.5 


S3 % 


13.2 








w% 


33.7 


40.5 


39.2 
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Fig. 7. Numerical simulation of the Xi Kou landslide: debris accumulation at different times and computed instantaneous flow rates 
(arrows). For clarity only rates exceeding 0.1 * qmax have been plotted. The contours of the reported slide are superimposed to the 
graph for reference. 
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Fig. 8. Numerical simulation of the Xi Kou landslide; water 
fraction (a) and class 3 fraction (b) in the body of the slide at 
t = 80, together with flow rates (arrows). For clarity only rates 
exceeding 0.1 * Qmax have been plotted. 



in Figs. 7-9. These figures show the debris height and 
the flow velocity vectors in the cells at different times, 
plotted over the planimetry of the site. 

The first panel of Fig. 7 shows the situation before 
the landslide. The debris is mainly located in the upper 
part of the slope and covers the platform. The following 
panels of Fig. 7 report intermediate flow snapshots. The 
simulated movement develops within the boundaries of 
the real event, branching to the secondary direction and 
interesting the hill as reported. The computed flow ve- 
locities are consistent with the preferential flow paths 
that result from the local morphology. 

It may also be seen (Fig. 8) how, at intermediate 
times, the model reproduces a characteristic front surge 
shape, with higher concentrations of the coarser mate- 
rial found at the boundaries of the active region. 

The larger discrepancies between the simulation and 
the event reported by Jiang Yun et al. (1991) are en- 
countered at the slope toe and along the secondary ton- 
gue. We note that they occur where the lack of data 
has forced us to resort to conjectures: at the toe, the 
concave shape of the valley is generated by the inter- 
polation, and is due to the lack of surrounding isolines; 
for the hill, the choice of modelling it as a cone of finer 
materials may have been arbitrary. We think however 
that the model satisfactorily captures the features of 
the event. A better calibration would be possible with 
more accurate and extensive data about the granulo- 
metric composition of the flowing mass. A complete 
calibration of the model would also require a statistical 



Fig. 9. Numerical hydrograph of the debris flow entering cell 
(25, 60): flow rate/base area and angle with the y-axis versus time. 



exploration of the parameters space for the maximum 
likelihood of the simulation. 

In the present case, the parameter values have been 
obtained empirically while trying to reproduce as better 
as possible the shape and the extension of the area in- 
vaded by the debris flow and the event duration. These 
values are reported in Table 4. The simplest square lat- 
tice connectivity was used here. 



5 Statistical properties 

The debris flow we model is essentially a dissipative pro- 
cess. Even though our model does not account for an 
energy balance, the global dissipative character is appar- 
ent when evaluating the potential energy of the material 
distribution, which is 



Wit)^gSj2pihtMht)(zii) 



h{i,t) 



(27) 



As can be seen from Fig. 10, W{t) is asymptotically 
decreasing, with small and irregularly spaced bursts of 
activity. These can be interpreted as temporary fluctu- 
ations, caused by the advection effect, which can occa- 
sionally draw some material counter-slope. Changes in 
the slope of W{t) are connected with the development 
of the different phases of the landslide, and are related 
to the actual debris distribution over the topography. 



Table 4. Parameter values for the Xi Kou slide simulation 



Q-i 


Qo 


Qi 


02 


03 


a 


13 


At 


0.030 


0.030 


0.013 


0.020 


0.040 


0.01 


0.03 


0.02 
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Fig. 10. Potential energy of the debris layer in time. The inset 
shows a zoom-in of the final part of the curve. 



Following the line of investigation of the sandpile lit- 
erature, we looked for a distinctive nonlinear behavior 
of the system by examining some of its statistical prop- 
erties. In order to determine if the observed behavior 
could be described properly as intermittent, as happens 
for other nonlinear systems, we followed the procedure 
outlined by Provenzale et al. (1993). This requires a 
multifractal analysis of some signal that accounts for the 
global behavior of the system in time. We considered, 
as representative of the whole process, the cumulative 
height M{t) of the material that is moved at each time 
step, i.e., 

b 3 

^W ^ imb-i(*' ^' *) + (1 + ") II %(«' «' t)] ■ (28) 

i v—1 k—0 

We presume that this quantity may be some sort of anal- 
ogous of the avalanche size for the sandpile models. We 
decided to refer to time series generated by the Xi Kou 
simulation, although the procedure might have been ap- 
plied as well to a completely idealized event. The time 
series of M(i), plotted both in lin-lin and in log-log scale 
is shown in Fig. 11. Both a power law trend, which is 
the effect of the asymptotic behavior mentioned in Sect. 
2.2.2, and the presence of irregular bursts are apparent 
in this plots. 

We analyzed two sub-series consisting of 32768 points 
of this signal, starting from t = 1000 and from t = 
1844.64. We estimated that, at these times, the most of 
the mass movement had already taken part (compare, 
e.g., the last two panels of Fig. 7), and the residual pro- 
cess could be considered nearly stationary. This is not 
entirely true, since the statistical properties of the two 
series are not found to be identical, and since a back- 
ground trend is still present. A detailed analysis shows, 
however, that this trend docs not affect the following 
conclusions. 

We then considered the numerical derivative D{t) = 
[M{t + At) - M{t)] /At, i.e., a fluctuation signal. This 



Fig. 11. Cumulative movement of material across the lattice 
during time, M{t). Insets (a) and (b) show details of the curve; 
inset (c) is a log- log plot of the same. Y-axis units arc of h/time. 



numerical differentiation is sufficient for filtering out the 
small and much slower relaxation trend present in the 
two sub-series. 

Fourier analysis of D{t) shows a noisy spectrum (Fig. 
12), without the clean power law decay exhibited by 
simpler sandpile models for the distribution of the ava- 
lanche sizes. 

Figure 13 shows the cumulative histogram of the fluc- 
tuations, which evidently follow a non-gaussian statistic. 
The values of the first two moments of the distributions 
and some other fine details slightly differ for the two sub- 
series. However, both histograms have a shape which is 
much closer to that of exponential distributions, with 
some extra features like a secondary bump for negative 
deviations. 

In parallel, we created two surrogate time series by 
phase randomization of the amplitude spectra of the 
original ones. The probability distributions of their fluc- 
tuations (Fig. 14), this time strikingly gaussian, denote 
that the nontrivial statistic is reflected in a phase lock- 
ing of the Fourier components. One conclusion of this 
fact is that the instantaneous transfers across the lat- 



Fig. 12. Amplitude spectra of D{t). Dotted lines are a power 
law fit of the last decade of the spectrum, included for reference. 
Here and afterwards (a) and (b) refer to the two sub-scries. 
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Fig. 13. Histograms of the fluctuations of the two sub-scries of 
M{t). The dotted parabolas correspond to gaussian distributions 
with the same first two moments. 



Fig. 15. Spectra of dq for the sub-series of M(t) (□) and for 
surrogates (o). Error bars represent 95% confidence intervals on 
the line fit. 



tice cannot be uncorrelated, and the landslide process 
has memory of its history. 

To characterize the intermittency we then computed 
the generalized fractal dimension dq of the sub-series, 
along with that of their surrogates. Following Proven- 
zale et al. (1993) and the references cited therein, we 
first computed the partition function 



Bs{e,q)^J2U 




(29) 



In the latter equation, d{t) is an appropriate measure 
function of the activity of the signal; numerical inspec- 
tion proved D{t)^ to be a suitable one. If the signal is 
multifractal, then, for small e, Bs{e,q) scales like 



Bs(e,q) ex e 



(q-l)d. 



(30) 



The generalized dimension dq is thus found for every q 
by a linear fit of Bs(e,q) in log- log coordinates, within 
the scaling range. The computation shows a good power 
law scaling of Bs{e,q) on the range At and T/30, that 
is over more than 3 decades. The multifractal spectra 
dq in the range < g < 4 are shown in Fig. 15. The 
surrogates have higher values of dq than those of the 
original series, and also a slower decrease of dq with 
increasing q. We therefore conclude that the original 
series have a self-affine behavior, which is peculiar since 
it is lost after the phase randomization; our signals M(t) 
can be then legitimately considered as multifractal. 



Fig. 14. Histograms of the fluctuations of the surrogate signals. 
The format is the same as in Fig. 13. 



6 Conclusions and perspectives 

We have set up a model of debris flow which is schematic 
and therefore computationally light. A preliminary vali- 
dation, based on the field data that was presently avail- 
able to us, showed that the model is able to capture 
most of the features of the natural process. In light 
of this agreement, it would be interesting to perfection 
the model and to calibrate it with more complete data, 
thus establishing a correspondence between its parame- 
ters and field measurements of geotechnical significance. 
A complete calibration shall base itself on a fully moni- 
tored and accessible landslide event, in which the granu- 
lometric composition at various locations and the extent 
of the slide are exactly known. A refinement of the flow 
rules may also be required. 

We have seen that the model behaves globally as a 
nonlinear system, giving rise to intermittent signals. In 
this our model is probably just an elaboration of the 
basic sandpile models, and follows their statistical char- 
acteristics. 

It would be interesting, but it falls beyond the scope 
of this paper, to make an explicit connection between 
our model and other lattice methods. The construction 
of our model was indeed suggested by the extensive ex- 
amples of lattice gases, but many important differences 
prevent us from making use, at this stage, of known re- 
sults for the latter. Some analogies can be drawn, and 
we just want to point at them, without attempting to 
establish the complete reasoning that connects lattice 
gases with continuum equations (Frisch et al., 1986). 

In lattice gases, self-excluding particles travel along 
the lattice branches, and collision rules that conserve 
mass, momentum and energy are implemented. Such 
particles may have different masses, or be of different 
species. In our model we have different moving ma- 
terials, but only material mass is conserved. The role 
of critical directions is perhaps equated to that of the 
gas molecules, in the sense that momentum is extracted 
from the cell in those directions. Our model may be per- 
haps compared in this context with lattice Boltzmann 
automata, because it employs real state variables rather 
than discrete ones. This is also suggested by the fact 
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that our automaton relies on a sort of local equilibrium 
assumption, since models a non-stationary process by 
assuming the local validity of an equilibrium flow law. 
Likewise, Boltzmann automata use equilibrium proba- 
bility distribution functions at each lattice site, and can 
already work with a single relaxation time form of the 
collision operator (Chen et al., 1992). 

Another parallel that could be made would be view- 
ing our landslide automaton as a sort of percolation 
model (see, e.g., Sahimi, 1993, for a review), in which 
the percolation network is continuously modified by the 
advected material. 

In our primary point of view, however, just a descrip- 
tive model was sought. In this aim, a finer resolution of 
the studied domain is expected to lead just to a greater 
accuracy of the model. In a lattice gas, instead, a finer 
resolution is considered in order to perform macro av- 
erages, which allow to recover continuum equations in 
the asymptotic limit. That might also be done with our 
set of cell rules, which are more complicated than the 
ones of the Navier-Stokes automata, perhaps in the di- 
rection of recovering some non-Newtonian debris flow 
equation. However, we note that even for the empiric 
constitutive relations there is no complete agreement in 
the literature, where generally just particular regimes of 
equilibrium flows are considered. In the general case, 
for flows with a free surface, varying slope, flow depth 
and flow rate, vertically averaged De Saint- Venant type 
equations are usually found. At the present state, we 
would not know to which continuum equation the model 
should tend in a thermodynamic limit. Moreover, our 
interest was in the large scale modelling of a landslide 
site, rather than on the fundamental form of the consti- 
tutive debris flow equations. 

Once better founded, we believe that this model could 
be not only a tool for a better simulation of hazardous 
events, but also open the possibility of numerical inves- 
tigation on more abstract topics, such as the connection 
between the statistical properties of the topography of 
arbitrary landscapes and that of the landslides taking 
part on them, or the response of unstable debris accu- 
mulations to random perturbations. 
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